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I.  INTRODUCTION 

To  gain  information  concerning  the  transient  response  of  a  complex 
structure  under  an  arbitrary  loading,  an  analysis  of  vibration  response  time 
history  data  is  required.  Unfortunately,  a  continuous  system  has  an  infinite 
number  of  degrees  of  freedom  from  which  data  can  be  extracted.  Taking 
measurements  at  all  of  possible  locations  on  the  structure  is  obviously  not 
feasible  due  to  constraints  on  time,  resources,  and  money.  Because  of  these 
constraints,  accelerometers  are  placed  at  a  few  specific  locations  on  a  structure 
and  the  data  is  then  analyzed.  While  this  information  is  quite  useful,  the 
collection  of  data  still  does  not  provide  a  clear  picture  of  the  response  at  the 
many  locations  where  data  was  not  taken.  More  importantly,  since  each  piece 
of  data  is  for  a  specific  location,  it  takes  quite  a  bit  of  imagination  to  get  an 
understanding  of  how  the  structure  as  a  whole  is  responding  to  the  loading. 

Due  to  the  limits  on  the  number  of  accelerometers  available,  the  shock 
trial  conducted  on  the  DDG-51  class  destroyer  provided  a  spatially  incomplete 
set  of  test  data.  In  addition  to  the  information  collected  during  the  test, 
knowledge  of  the  transient  response  at  locations  on  the  ship  where  data  was 
not  taken  was  desired.  By  expanding  this  collection  of  test  data  into  a  spatially 
complete  set,  the  transient  response  could  be  visualized  in  a  computer 
simulation. 

This  thesis  investigates  the  strengths  and  weaknesses  of  expanding 
spatially  incomplete  test  data  using  finite  element  model  reduction  methods. 
These  reduction  methods  are  used  to  create  a  transformation  matrix  which 
makes  it  possible  to  expand  a  relatively  small  collection  of  simulated  test  data 
into  a  set  that  corresponds  to  a  much  larger  number  of  degrees  of  freedom. 


The  transformation  matrix  is  created  from  a  finite  element  model  of  the 
structure  whose  system  mass  and  stiffness  matrices  have  been  partitioned  in 
accordance  with  the  locations  of  the  accelerometers  on  the  actual  structure. 
The  final  collection  of  data  becomes  the  time  histories  of  all  the  coordinates 
in  the  full-order  finite  element  model  and  has  the  same  number  of  degrees  of 
freedom. 

The  full-order  set  of  data  is  then  animated  to  create  a  real-time 
visualization  of  the  dynamic  response  of  the  structure  as  a  continuous 
system.  Every  column  in  the  expanded  matrix  is  a  representation  of  the 
dynamics  of  each  node  in  the  finite  element  model  at  a  specific  time  interval. 
With  this  information,  a  very  clear  understanding  of  the  vibration  response 
of  a  complex  structure  is  possible. 

There  are  two  different  reduction  methods  used  in  this  thesis  to  create 
the  transformation  matrices.  The  first  is  the  static,  or  Guy  an,  reduction  and 
the  second  is  the  Improved  Reduced  System  (IRS)  reduction.  The  IRS 
method,  while  initially  more  computationally  expensive,  provides  an 
improvement  because  it  approximates  the  inertia  forces  that  are  present  in  a 
dynamic  problem.  A  sensitivity  analysis  for  each  simulation  is  provided  to 
highlight  the  accuracy  of  both  the  static  and  the  dynamic  reduction  methods. 


II.  THEORY 

A.         MATRIX  PARTITIONING 

In  the  physical  coordinate  system,  the  expansion  of  the  collected  time 
history  data  can  be  expressed  in  the  following  way: 


M  (D 


where  {xaj  is  a  vector  containing  the  dynamic  response  data  collected  from 
the  structure,  {x0}  is  the  response  at  coordinates  other  than  the  data  collection 

points,  and  [T]  is  a  non-square  transformation  matrix  used  in  the  expansion 
of  {xa}  to  the  full-order  response  set,  and  [I]  is  the  identity  matrix. 

The  subscript  'a'  refers  to  the  analysis  set  of  model  coordinates,  which 
from  here  on  will  be  called  simply  the  'a-set'.  This  is  the  set  of  coordinates  in 
the  model  that  corresponds  to  the  locations  on  the  structure  where  data  was 
taken.  Conversely,  the  subscript  'o'  refers  to  the  omitted  set.  This  'o-set'  refers 
to  all  subsequent  model  coordinates  where  data  is  not  taken.  By  operating  on 
the  'a-set'  time  history  data  set  with  the  proper  transformation  matrix,  the  'o- 
set'  dynamic  response  is  calculated. 

Generally,  due  to  the  size  of  complex  models,  the  a-set  system  will  be 
much  smaller  than  the  o-set  system.  The  model  of  a  complex  structure  will 
normally  contain  a  huge  number  of  degrees  of  freedom  while  the  number  of 
accelerometers  is  relatively  small  in  comparison. 
B.         STATIC  TRANSFORMATION  MATRIX 

By  partitioning  the  stiffness  relation,  f  =  kx,  the  static  transformation 
matrix  can  easily  be  found.  The  equation  becomes: 


K„      K 


(2) 


By  taking  care  to  partition  the  equation  such  that  there  are  no  excitations 
applied  at  the  o-set  coordinates,  then  {fn}  can  be  set  equal  to  {0}.  Now  the 

relation  between  the  o-set  and  a-set  response  can  be  expressed  as: 


k}=[-KX.]k} 


(3) 


The  static  transformation  matrix  is  then  [Ref.  1], 


T     = 

stat 


I 

-k:'k 


(4) 


C         IRS  TRANSFORMATION  MATRIX 

The  derivation  of  the  IRS  reduction  transformation  matrix  begins  with 
the  partitioned  equation  of  motion  for  a  vibrating  system  [Ref.  2]: 


ML    M, 
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x. 


(5) 


Assuming  a  harmonic  motion  of  frequency,  Q,  the  acceleration  term  in 
equation  (5)  can  be  replaced  with  x  =  -Q.2x .   Again  partitioning  such  that  there 
are  no  excitations  at  the  o-set  coordinates,  we  arrive  at  the  exact  structural 
dynamics  reduction  relation  [Ref.  3]: 


{x„}  =  [l-^K:X0]-'[-K:Xa  +Q2K0Moa]{xa}  (6) 

(a)  (b) 

Ideally,  the  desired  transformation  matrix  will  be  frequency  independent  and 
Eq.  (6)  is  dependent  on  the  driving  frequency  by  the  Q.2  term.  Truncating  the 
binomial  expansion  of  Eq.  (6a)  after  the  terms  that  include  Q2  yields: 

{x0}[-k;X  +^2(k;X.  -k;XXXJ]W  (?) 

The  frequency  dependency  is  then  eliminated  by  the  following  approximation 
of  the  acceleration: 

^2{xfl}  =  [M,fJ-I[KrtJ  (8) 

where  the  statically  reduced  mass  and  stiffness  matrices  are  given  by  the 
following  [Ref.  2,3]: 

MM  =  TlMT_  (9a) 

K^TlKT^,  (9b) 

Substitution  of  Eq.  (8)  into  Eq.  (7)  provides: 

K}[-k;Xu  +[k;X„  -  k;X^k  Jm^-^k  J]{xj         (io) 

which  is  the  IRS  reduction  relationship  between  the  a-set  and  o-set 


coordinates.  The  IRS  reduction  transformation  matrix  can  then  be  written  as 
[Ref.  2,3]: 


T     = 


-k;X„  +[k;>ou  -k;:mook;:kou]m;:k, 


(ii) 


D.         VISUALIZATION  MATRIX 

After  operating  on  {xa}  with  the  transformation  matrix,  the  response 

for  the  entire  structure  is  partitioned  in  the  form: 


{x(0}  = 


'{xm}   .   .   {xa(rj}' 
{%M}   .   .   k(rj} 


(12) 


where  the  columns  in  the  matrix  are  a  representation  of  the  dynamic 
response  at  each  node  in  the  model  at  a  specific  instant  in  time.  However,  this 
partitioning  does  not  correspond  to  the  nodes  in  the  model,  therefore  Eq.  (12) 
must  again  be  partitioned  such  that  it  is  in  line  with  the  nodes  in  the  model. 
To  get  a  real-time  visualization  of  the  structure,  simply  march  through  the 
matrix  taking  snapshots  of  the  model  at  each  time  interval. 


III.  FINITE  ELEMENT  MODEL  FORMULATION 

The  finite  element  method  used  in  the  creation  of  the  plate  model  for 
this  thesis  is  based  on  the  shear  deformable  displacement  formulation.  The 
element  chosen  is  the  four-noded  quadrilateral  plate  element.  Each  node  in 
the  model  has  three  degrees  of  freedom  which  are  the  transverse 
displacement,  w  ,  and  the  two  rotations  about  the  midplane,  6X  and  6y.  Using 

the  following  bilinear  isoparametric  shape  functions  [Ref.  4], 


»,(&!?)  =  j0  -m-V)  d3a) 

H2{S,n)  =  1(1 +  £)(!- 77)  (13b) 

ff,({,ij)  =  1(1  +  9(1 +  77)  (13c) 

H4fari)  =  1(1  -«)(!  +  77)  (13d) 


the  transverse  displacement  and  slopes  are  interpolated  as: 

w  =  £/7,(£77)w,  (14a) 

(=1 

ex  =  tH^Ml  (14b) 

1=1 

sy  =  lm-n)(eyl  (ho 


;=i 


These  isoparametric  shape  functions  are  defined  in  terms  of  a  normalized 

domain  -1  <£<1  and  -1  <r\<  1. 

A.        INTEGRATION  TECHNIQUE 

As  previously  stated,  the  shape  functions  used  in  this  model  are  for  the 


bilinear  isoparametric  element.  This  type  of  element  was  chosen  in  order  to 
make  the  finite  element  model  general  in  the  sense  that  it  can  be  applied  to 
any  plate  geometry.  The  elements  in  the  physical  coordinate  system  are 
mapped  into  the  isoparametric  coordinate  system  by  the  following  relations: 


(15a) 
(15b) 


1=1 


where  x  and  y  are  physical  coordinates  corresponding  to  the  nodes  in  the 
element.  Gauss-Legendre  quadrature  is  used  to  perform  all  integrations. 
B.         ELEMENT  MATRIX 

1.  Elemental  Stiffness  Matrix 

Without  providing  the  derivation,  the  elemental  stiffness  matrix,  [K*], 

for  plate  bending  is  [Ref.  4]: 


[K<]  =  ^  JBTbJ)bBbdQ  +  OiJBlDfi/iQ 


a' 


cic 


(16) 


where  h  is  the  thickness  of  the  plate,  £  equals  —  and  is  the  shear  energy 

6 

correction  factor,  Qe  is  the  two  dimensional  element  domain,  and 
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(18) 

The  constitutive  equation  for  the  plane  stress  condition  is: 


[D>]=ii 


1  v 
v  1 
0    0 


0 

0 
1-v 


(19) 


while  the  constitutive  equation  for  shear  is, 


[D,]  = 


G     0 
0     G 


(20) 


G  is  the  shear  modulus,  E  is  the  modulus  of  elasticity,  and  v  is  Poisson's 
ratio.  One  thing  of  importance  to  be  noted  about  Eq.  15  is  that  the  shear  energy 
becomes  dominant  over  the  bending  energy  as  the  thickness  of  the  plate 
becomes  small  relative  to  the  side  length.  This  is  called  shear-locking.  To 
account  for  this  problem,  a  reduced  integration  technique  is  used.  The 
bending  term  is  integrated  exactly  using  2x2  Gauss-Legendre  quadrature  while 
the  shear  term  is  under-integrated  using  lxl  Gauss-Legendre  quadrature. 

2.  Elemental  Mass  Matrix 

The  elemental  mass  matrix  is  found  by  simply  integrating  the  properly 
weighted  shape  functions.  The  general  equation  for  determining  this  matrix 


is: 


[W]=\PA 


a' 


[Hx     H2     H2     H4]dQ. 


(21) 


In  Eq.  (19),  P  is  the  density  of  the  plate  element,  A  is  the  elemental  area,  and 
[H]  is  the  same  shape  functions  as  in  Eq.  (12).  2x2  Gauss-Legendre  quadrature 
is  used  to  perform  all  integrations  for  the  element  mass  matrix. 
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IV.  COMPUTER  SIMULATIONS 

The  investigation  of  expanding  time  history  data  was  conducted  with 
the  aid  of  computer  simulations.  The  dynamic  response  at  each  node  of  a 
plate  finite  element  model  was  used  to  simulate  the  'actual'  solution.  From 
this  baseline  model,  an  a-set  system  of  time  histories  was  chosen  which  is 
assumed  to  be  the  data  taken  from  an  actual  structure.  The  system  matrices  in 
the  finite  element  model  were  then  partitioned  in  order  to  create  the 
transformation  matrix  required  to  expand  the  a-set  time  histories  to  full- 
order.  By  comparing  the  dynamic  response  of  the  finite  element  model  to  the 
response  found  by  expanding  the  a-set  system,  the  validity  of  the  process 
could  be  studied  and  verified. 
A.        EXAMPLE  (1):  TWO  DOF  IN  THE  A-SET 

The  initial  simulations  are  conducted  on  a  square,  flat  plate  (see  Figure 
1)  comprised  of  twenty-five  elements  and  thirty-six  nodes.  Degrees  of  freedom 
42  and  87  are  assumed  to  be  the  'actual'  response  and  therefore  make  up  the  a- 
set  system.  Conversely,  the  o-set  system  is  comprised  of  all  other  DOF.  The 
two  DOF  in  the  a-set  system  correspond  to  the  transverse  motion  at  nodes  14 
and  29.  The  external  force  is  a  blast  function  (see  Figure  2)  applied  at  node 
twenty-one.  Each  external  spring  has  a  stiffness  equal  to  1000  lb-in  and  is 
attached  to  the  ground. 

Two  separate  simulations  are  performed  using  this  plate  and  a-set 
system.  The  first  simulation  uses  the  static  transformation  matrix  to  expand 
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the  a-set  and  the  second  uses  the  IRS.  Both  are  then  compared  to  the  full- 
order  plate  response.  Unfortunately,  it  is  impossible  to  show  the  visualization 
on  paper  so  the  comparison  will  be  on  a  per  DOF  basis. 
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Figure  1  Plate  Model 
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Figure  2  Forcing  Function 

1.  Full-Order  Plate  Results 

The  dynamic  response  of  the  full-order  plate  model  was  solved  by- 
transforming  the  general  equation  of  motion, 


[M]{x}  +  [K]{x}  =  {f(t)} 


(22) 


into  modal  coordinates  by  the  following  relation: 


13 


Wt)}  =  fo]{q(t)}  (23) 

in  which  {q(t)}  is  the  modal  displacements  (as  a  function  of  time)  and  [(j>]  is  a 
matrix  of  mass  normalized  modeshapes  found  from  solving  the  typical 
structural  dynamics  eigenvalue  problem: 

[K-^2M]W  =  {0}  (24) 

The  modal  displacements  are  found  by  numerically  solving  a  convolution 

integral  [Ref.  5].  The  damping  in  the  plate  is  assumed  to  be  2%  and  is  inserted 

into  the  integral  calculation.  Once  {q  (t)}  is  found,  Eq.  (21)  is  then  used  to 

convert  back  to  the  physical  coordinate  system. 

Now,  the  time  history  at  every  node  in  the  plate  is  known.  These  results  will 

be  used  to  compare  the  accuracy  of  the  static  and  dynamic  expansions  of  the  a- 

set. 

2.  Static  Transformation  Results 

Figures  III,  IV,  and  V  are  per-node  samples  of  the  full-order  plate 
model  response  in  comparison  to  the  response  found  by  the  static  expansion 
of  the  2  DOF  a-set.  It  is  clear  that  the  static  transformation  will  provide  an 
adequate  approximation  of  the  'true'  response  at  certain  nodes. 
Unfortunately,  Figure  3  proves  that  the  difference  between  the  'true'  solution 
and  the  expanded  set  can  be  greater  than  50%. 
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Figure  3  Static  Response  at  Node  3 
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Figure  4  Static  Response  at  Node  17 
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Figure  5  Static  Response  at  Node  27 


Having  errors  of  this  magnitude  can  lead  to  very  misleading  results  in 
the  final  visualization.  It  will  provide  an  idea  of  how  the  structure  vibrates 
but  the  animation  is  distorted  and  skewed.  If  a  more  accurate  representation 
is  the  requirement,  then  other  steps  must  be  taken.  Obviously,  one  way  to  get 
a  more  accurate  solution  is  to  increase  the  size  of  the  a-set  system.  An 
example  solution  of  the  larger  a-set  is  provided  later.  The  drawback  to  this  is 
more  data  must  be  taken  which  is  oftentimes  impossible  due  to  limits  on  the 
number  of  accelerometers  available.  Another  way  to  change  the  accuracy  of 
the  solution  is  to  alter  the  locations  where  the  data  is  actually  taken  on  the 
structure.  An  example  is  also  provided  later  in  this  chapter  of  greatly 
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decreasing  the  accuracy  of  the  static  solution  by  choosing  an  a-set  system 
badly. 

a.  Error  Analysis 

To  understand  the  differences  in  the  full-order  and  static  plate 
solutions  without  providing  the  time  histories  at  every  node,  a  sensitivity 
analysis  was  conducted. 

Sensitivity  Analysis  Between  Static  and  Full-Order  Response 
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0     0 


Figure  6  Error  of  Static  Reduction  Method 
note:  Figure  is  orientated  such  that  the  0/0  point  is  node  1 


The  largest  difference  between  the  full-order  plate  response  and 

the  expanded  time  history  set  was  found  and  plotted.  This  analysis  is 

essentially  a  'worst-case'  scenario  because  the  code  searched  through  time  to 
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find  the  absolute  maximum  error  at  each  node  and  these  errors  did  not 
always  occur  at  the  same  instant  in  time. 

The  maximum  differential  appears  to  be  approximately  0.1 
inches,  but  the  movement  of  the  plate  is  only  about  0.3  inches.  Under  these 
circumstances,  an  error  of  this  magnitude  is  quite  significant,  and  therefore 
the  static  reduction  probably  should  not  be  used. 

3.  IRS  Transformation  Results 

Using  the  existing  plate  model  and  forcing  function  (see  Figures  1  and 
2),  the  solution  for  the  IRS  expanded  data  is  provided  as  before.  The  same 
nodes  are  given  as  samples  here: 
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Figure  8  IRS  Response  at  Node  17 
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Figure  9  IRS  Response  at  Node  27 
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These  figures  show  that  the  IRS  expanded  data  is  much  more  accurate 
in  comparison  to  the  static  transformation.  Even  with  an  a-set  of  only  2  DOF 
out  of  the  possible  108,  the  IRS  transformation  provides  a  very  good  solution. 
Clearly,  the  correction  in  the  IRS  transformation  matrix  due  to  inertial  forces 
will  improve  the  accuracy  of  the  solution  immensely. 

a.         Error  analysis 

A  sensitivity  analysis  for  the  IRS  expanded  data  is  also  provided: 


Sensitivity  Analysis  Between  IRS  and  Full-Order  Response 
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Figure  10  Error  of  IRS  Reduction  Method 
note:  Figure  is  orientated  such  that  the  0/0  point  is  node  1 

The  maximum  difference  in  Figure  10  is  less  than  0.02  inches. 
This  accuracy  is  an  order  of  magnitude  better  than  the  static  reduction.  The 
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approximation  of  the  inertial  forces  inherent  in  the  IRS  transformation 

matrix  appears  to  work  quite  well. 

B.         EXAMPLE  (2):  FOUR  DOF  IN  THE  A-SET 

This  example  shows  that  the  accuracy  in  both  solutions  can  be  greatly 
improved  by  increasing  the  size  of  the  a-set.  By  merely  increasing  the  size  of 
the  a-set  by  two,  it  will  be  seen  that  the  improvement  in  accuracy  is  at  least  an 
order  of  magnitude.  As  mentioned  earlier,  although  desirable,  increasing  the 
number  of  data  points  may  not  be  feasible. 

For  this  simulation,  the  forcing  function  used  is  the  same  as  in 
Figure  2.  The  plate  model,  shown  in  Figure  1,  is  changed  such  that  the  data  is 
instead  extracted  at  nodes  8, 11,  26,  and  29.  In  the  transverse  direction,  these 
nodes  correspond  to  degree  of  freedom  24,  33,  78,  and  87.  All  other  variables 
in  the  simulation  are  unchanged. 

On  a  per  node  basis,  the  results  are  very  hard  to  comprehend  because 
the  graphs  are  almost  exactly  on  top  of  each  other.  For  the  IRS  expanded  data, 
no  difference  can  be  detected.  Because  of  this,  only  the  sensitivity  analysis  is 
provided. 

1.         IRS  and  Static  Transformation  Error  Analysis 

By  viewing  these  graphs,  it  becomes  clear  how  the  accuracy  can  by 
greatly  increased  by  adjusting  the  number  of  DOF  in  the  a-set  system.  The 
greatest  difference  in  inches  of  the  response  for  the  static  solution  is  0.01. 
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Sensitivity  Analysis  Between  Static  and  Full-Order  Response 
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Figure  11  Error  of  Static  Reduction  Method 
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Figure  12  Error  of  IRS  Reduction  Method 
note:  Figures  are  orientated  such  that  the  0/0  point  is  node  1 
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The  difference  for  the  IRS  solution  is  much  smaller  at  0.00015  in.  As  expected, 
the  IRS  reduction  method  remains  the  more  accurate  of  the  two,  but  both 
methods  provide  very  good  solutions.  If  the  visualization  could  be  shown 
here,  one  would  not  be  able  to  discern  a  difference  between  the  two 
animations  and  the  'true'  solution. 
C         EXAMPLE  (3):  A-SET  INCLUDES  CONSTRAINED  NODES 

For  this  simulation,  the  accelerometers  are  assumed  to  be  located  at 
nodes  which  are  constrained  to  ground  through  external  springs.  For  the 
plate,  these  nodes  are  numbered  1,  6,  31,  and  36.  The  corresponding  DOF's  in 
the  a-set  system  are  3, 18,  93,  and  108,  respectively. 

1.         Static  Transformation  Results 

Figures  13  and  14  are  a  sampling  of  the  results  from  this  particular 
simulation  and  they  show  a  very  interesting  situation  to  consider.  These  time 
history  samples  clearly  indicate  that  the  static  transformation  matrix  will  not 
provide  an  accurate  solution  for  this  particular  a-set.  On  closer  inspection  of 
the  two  previous  graphs,  it  is  seen  that  the  static  transformation  provided  the 
exact  same  solution  at  nodes  10  and  29.  In  fact,  although  not  shown,  this  is  the 
exact  response  calculated  at  every  node  in  the  plate,  and,  when  placed  in 
animation,  the  plate  exhibits  rigid  body  motion  in  the  transverse  direction. 

The  reason  for  this  rigid  body  motion  can  be  attributed  to  the  location 
of  the  a-set  system.  By  taking  time  history  data  at  only  those  nodes 
constrained  to  ground  through  springs,  the  transformation  matrix  imparts 
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Figure  13  Static  Response  at  Node  10 
Time  History  at  Node  29 
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Figure  14  Static  Response  at  Node  29 
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the  strain  energy  of  the  springs  to  the  system  and  ignores  any  of  the  plate's 
internal  strain  energy.  In  other  words,  the  response  calculated  in  the  o-set  is 
effectively  'zeroed  out'  by  the  transformation  matrix.  By  discarding  the 
motion  of  the  individual  o-set  nodes,  the  plate  moves  with  the  springs  as  a 
rigid  body.  A  detailed  explanation  of  this  behavior  is  provided  in  the  next 
chapter. 

a.  Error  analysis 

The  inaccuracy  of  the  solution  becomes  clear  with  the  sensitivity 
analysis. 


Sensitivity  Analysis  Between  Static  and  Full-Order  Response 
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Figure  15  Error  of  Static  Reduction  Method 
note:  Figure  is  orientated  such  that  the  0/0  point  is  node  1 
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Obviously,  rigid  body  motion  is  an  inaccurate  solution  to  this 
situation.  As  expected,  the  nodes  at  the  four  corners  exhibit  no  error  while  the 
center  of  the  plate  has  the  greatest.  From  Figure  13,  the  maximum  amplitude 
of  the  plate  motion  is  about  0.3  in.  which  is  the  same  as  the  maximum  error 
shown  here.  Example  (2)  proved  that  a  very  accurate  solution  can  be  found 
using  the  static  reduction  method  when  the  a-set  contains  four  DOF.  This 
example  shows  that  care  must  be  taken  when  using  the  static  reduction 
method  that  the  accelerometers  are  not  placed  at  constrained  nodes. 

2.  IRS  Transformation  Results 

The  problems  using  the  static  reduction  method  are  not  exhibited  by 
the  IRS  method.   Due  to  the  correction  in  the  derivation  of  the 
transformation  matrix  due  to  the  inertial  forces,  the  motion  of  the  plate  is 
taken  into  account.   For  this  reason,  the  actual  response  of  the  plate  is 
calculated. 

The  following  time  history  samples  prove  the  usefulness  of  the  IRS 
reduction  method  when  data  must  be  taken  at  nodes  constrained  to  ground 
through  a  spring: 
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Figure  16  IRS  Response  at  Node  10 
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a.         Error  analysis 

Figure  18  has  the  same  basic  shape  as  Figure  15,  but  the 
maximum  error  this  time  is  only  0.02  in.  The  inertial  force  correction  to  the 
IRS  transformation  matrix  enables  a  very  accurate  solution  to  be  found  even 
though  the  strain  energy  of  the  plate  is  ignored.  For  this  reason,  using  the  IRS 
method  allows  for  fewer  limitations  in  the  location  of  the  a-set  system. 


Sensitivity  Analysis  Between  IRS  and  Full-Order  Response 
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Figure  18  Error  of  IRS  Reduction  Method 
note:  Figures  are  orientated  such  that  the  0/0  point  is  node  1 
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V.  LESSONS  LEARNED 

The  three  examples  provided  in  the  previous  chapter  are  only  a  few  of 
the  simulations  conducted  in  the  investigation  of  expanding  spatially 
incomplete  time  histories  to  full-order.  Several  other  situations  requiring 
investigation  were  checked  and  will  be  addressed  here.  This  chapter 
highlights  the  lesson's  learned  from  all  the  computer  simulations  without 
providing  the  detail  seen  previously. 

A.  ACCURACY  DEPENDS  UPON  THE  NUMBER  OF  DOF  IN  THE  A-SET 
As  expected,  the  number  of  DOF  in  the  a-set  system  greatly  affects  the 

accuracy  of  the  solution.  The  increase  in  accuracy  seen  in  changing  the  a-set 
system  from  Example  (1)  to  Example  (2)  is  excellent  proof  of  this.  Actually,  a 
dramatic  improvement  could  be  seen  by  increasing  the  number  of  DOF  in  the 
a-set  to  three.  Interestingly,  independent  of  the  size,  shape,  or  number  of  DOF 
in  the  model,  a  system  of  only  one  DOF  in  the  a-set  did  not  provide  a  viable 
solution  for  either  method.  Therefore  the  minimum  number  of  DOF  in  the 
a-set  is  restricted  to  two.  Simulations  were  not  run  for  situations  where  the 
number  is  greater  than  four  because,  as  seen  in  Example  (2),  an  a-set 
containing  four  DOF  expands  into  a  solution  that  is  essentially  exact. 

B.  LOCATION  OF  THE  A-SET  IS  CRITICAL 

In  addition  to  the  number  of  DOF,  the  accuracy  of  the  solution  can  be 
greatly  affected  by  the  location  of  the  a-set.  Two  general  cases  have  been  found 
that  should  be  taken  into  account.  The  first  is  to  pick  an  a-set  location  that  is 
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certain  to  impart  strain  energy  into  the  transformation  matrix  in  order  to 
avoid  the  rigid-body  modes.  This  restriction  is  of  vital  importance  if  the  static 
reduction  is  chosen.  The  second  is  to  ensure  that  the  time  history  data  is 
taken  at  points  'distant'  from  one  another. 

Example  (3)  revealed  a  problem  with  the  static  reduction  method  when 
the  nodes  chosen  for  the  a-set  are  constrained.  As  mentioned  earlier,  the 
reason  is  that  the  only  strain  energy  imparted  into  the  system  is  the  stiffness 
of  the  springs  while  the  internal  strain  energy  of  the  plate  is  'zeroed  out'. 
Mathematically,  this  idea  can  be  represented  by  first  partitioning  the  stiffness 
matrix  in  the  following  manner: 


[K]  = 


KL+KL 


■"■^•oo 

Kp 

*     00 


(25) 


where  the  superscript  p  is  used  to  indicate  the  plate  and  s  refers  to  the 
springs.  Let  x  be  a  displacement  vector  that  imparts  no  internal  strain  energy. 
This  now  defines  the  a-set  and  the  o-set.  Further  partitioning  of  Eq.  (23)  yields: 


[K]= 


Kp 

aa 

Kp 


+ 


K 
o 


(26) 


Now,  applying  Eq.  (24)  to  the  general  stiffness  relation,  f  =  kx,  the  following 
relation  between  the  a-set  and  the  o-set  is  derived  much  in  the  same  manner 
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as  Eq.  (3): 

M  =  -KVK]{xa}  (27) 

By  condensing  the  stiffness  of  the  grounded  springs  out  of  the  static 
transformation  equation,  the  o-set  exhibits  the  response  of  an  unrestrained 
structure.  The  rigid  body  motion  exhibited  in  Example  (3)  is  a  result  of  the 
stiffness  matrix  being  partitioned  in  this  manner. 

A  problem  that  affects  the  accuracy  of  both  the  static  and  IRS  reduction 
methods  is  choosing  the  a-set  system  at  points  that  are  'too  close'  together.  In 
this  situation,  common  sense  must  be  depended  upon  because  the  concept  of 
'close'  and  'far'  depends  on  the  size  and  geometry  of  the  structure.  Several 
simulations  were  conducted  to  study  this  and  the  best  solution  is  simply  to 
ensure  the  accelerometers  are  spread  out  to  locations  that  best  cover  the  entire 
structure.  Having  said  this,  it  is  important  that  the  a-set  include  locations 
where  the  motion  of  the  structure  has  the  greatest  amplitude.  In  the  case  of 
the  plate  model,  the  amplitude  of  the  vibration  is  greatest  at  the  center  of  the 
plate  so  the  a-set  should  include  time-history  information  from  this  area.  A 
simulation  was  run  where  the  a-set  was  comprised  of  edge  nodes.  The  plate 
moves  very  little  at  its  edges  and  subsequently  the  solution  was  inaccurate  for 
both  methods.  Sound  engineering  judgment  must  be  employed  when 
choosing  an  a-set  because  it  is  critical  to  getting  an  accurate  solution. 
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C         ROTATIONS  AS  PART  OF  THE  A-SET 

In  general,  using  rotations,  or  a  combination  of  rotations  and 
translations  as  part  of  the  a-set  does  not  seem  to  alter  the  accuracy  of  the 
solution.  The  same  rules  apply  for  rotations  as  for  the  transverse 
displacements.  A  larger  a-set  still  returns  a  more  accurate  solution  and  the 
location  is  still  important.  In  this  situation,  however,  being  closer  to  the  edge 
produced  more  accurate  solutions.  This  makes  sense  because  the  edges  are  the 
locations  of  the  largest  rotations.  Using  rotations  in  the  a-set  is  an  advantage 
if  applying  the  static  reduction  method  because  there  is  no  longer  a  restriction 
on  the  use  of  nodes  constrained  through  translational  springs. 

While  using  rotations  is  not  a  problem  computationally,  care  must  be 
taken  to  ensure  that  the  structure  is  rotating  along  the  axis  corresponding  to 
the  chosen  DOF.  If  not,  including  these  DOF  in  the  a-set  serves  no  purpose  as 
no  strain  energy  is  imparted  into  the  system  and  the  plate  will  again  exhibit 
rigid  body  motion. 
D.        VISUALIZING  A  STRUCTURE 

Animation  in  MATLAB  is  very  computer  memory  intensive.  A 

simple  structure  such  as  the  plate  can  tax  the  memory  of  most  computers  and 

cause  the  program  to  crash.  As  such,  a  very  large  structure,  or  one  with  many 

nodes,  would  be  almost  impossible  to  visualize.  To  aid  in  overcoming  this 

problem,  the  code  should  direct  the  computer  to  reduce  the  size  of  the 

graphics  window  for  the  animation.  Although  this  does  not  allow  for  as  nice 

a  visualization,  memory  is  saved  for  computation  instead  of  being  allocated 
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to  the  graphics  window.  By  doing  this,  the  visualization  can  be  of  longer 
duration  thereby  resulting  in  a  better  understanding  of  the  dynamics  of  the 
structure. 

A  better  tool  to  use  in  performing  a  dynamic  visualization  is  I-DEAS. 
Unfortunately,  certain  versions  will  not  allow  the  user  to  expand  time 
histories.  An  outline  can  be  found  in  the  appendix  that  describes  the  steps 
needed  to  perform  an  animation  in  I-DEAS. 


33 


34 


VI.  CONCLUSIONS  AND  RECOMMENDATIONS 

The  main  objective  of  this  thesis  was  to  investigate  the  process  of 
visualizing  transient  structural  response  by  expanding  spatially  incomplete 
time  history  data  to  full-order.  The  response  at  every  node  of  a  full-order 
plate  model  was  solved  which  provided  a  base  model  for  comparison.  From 
this  collection  of  nodal  responses,  a  few  of  the  time  histories  were  extracted 
which  simulated  the  experimental  data  taken  from  an  actual  structure.  After 
choosing  the  a-set,  the  finite  element  system  matrices  were  partitioned  to 
calculate  the  appropriate  transformation  matrix.  The  expanded  data  was  then 
compared  to  the  full-order  response  in  the  effort  of  evaluating  the  accuracy  of 
the  process. 

A.        CONCLUSIONS 
The  simulations  have  shown  the  following: 

1.  To  get  an  accurate  solution,  the  minimum  number  of  DOF's  allowed 

to  make  up  the  a-set  system  on  a  flat  plate  is  two.  Provided  the 
accelerometers  are  properly  placed,  an  a-set  comprised  of  four  DOF 
will  produce  a  very  accurate  response  for  both  the  static  and  the  IRS 
methods. 

2.  The  placement  of  the  accelerometers  is  very  important.  If  the  static 

reduction  method  is  chosen,  the  data  must  be  taken  at  nodes  certain 
to  produce  internal  strain  energy.  Regardless  of  method,  the  data 
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must  be  extracted  at  locations  'far'  enough  from  each  other  to  get  a 
sample  of  the  entire  structure,  and  the  data  collector  should  try  to 
pick  points  on  the  structure  where  vibration  is  greatest. 

3.  The  same  restrictions  that  apply  to  using  translations  in  the  a-set  can 
be  applied  to  rotations. 

4.  The  visualization  process  is  very  memory  intensive  therefore  clever 

programming  is  required  to  maximize  the  memory  available  for 
graphics. 

Visualizing  transient  response  provides  a  unique  insight  on  the 
motion  of  a  vibrating  structure.  A  much  better  understanding  of  the 
dynamics  is  gained  when  a  designer  can  visualize  the  motion  of  a  spatially 
complete  body  instead  of  the  time  histories  of  a  few  select  nodes. 
B.         RECOMMENDATIONS 

This  study  outlines  a  few  of  the  strengths  and  weaknesses  inherent  to 
the  Guyan  and  IRS  reduction  methods  in  the  attempt  to  animate  transient 
response.  While  many  simulations  were  explored,  the  investigation  was 
rather  limited  in  scope.  In  addition  to  the  Guyan  and  IRS  methods,  there  are 
many  reduction  methods  that  could  be  explored.  To  validate  the  process, 
experimental  data  must  be  loaded  into  the  program  and  studied.  A  similar 
study  should  be  directed  towards  visualizing  the  response  of  three 
dimensional  structures.  In  addition,  it  would  be  useful  to  investigate 
animating  transient  response  in  six  degrees  of  freedom. 
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The  accuracy  of  this  process  is  heavily  dependent  on  the  location  of  the 
a-set.  There  are  analytical  techniques  available  that  help  find  the  'best'  a-set 
but  ordinarily  these  are  used  when  applying  model  reduction  methods  to  find 
modeshapes.  These  techniques  should  to  be  investigated  to  prove  if  they  can 
also  be  applied  to  this  process. 

The  motivation  for  this  thesis  will  be  realized  when  the  time  history 
data  of  the  DDG-51  class  shock  trial  is  loaded  into  I-DEAS,  expanded  to  full- 
order,  and  then  animated  in  a  transient  response  visualization. 
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APPENDIX  A.  EXPANDING  TIME  HISTORIES  IN  I-DEAS 

Step-by-step  procedures  for  expanding  transient  time  history  data  using 
I-DEAS  software  is  provided.  This  option  is  not  available  using  I-DEAS 
Master  Series  3  if  the  software  is  set  up  to  run  on  Hewlett-Packer  computers. 
A.  ACCESS  I-DEAS  SIMULATION  INTERFACE  SOFTWARE 

1.  Master  Modeler  /Surfacing  Task  -  load /build  FE  model 

2.  Meshing  Task  -  mesh  model  and  enter  material  properties 

3.  Boundary  Condition  Task  - 

a.  Choose  normal  mode  dynamics 

b.  Apply  restraints 

c.  Create  boundary  condition  set/ Choose  normal  mode 
dynamics  -  Guyan 

d.  Choose  analysis  set  (a-set) 

4.  Model  Solution  Task 

a.  Create  solution  set 

b.  Initiate  solve 

c.  Access  I-DEAS  model  response  software 

5.  Model  Response  Software 

a.  Excitation  definition 

-  Create  force /displacement  excitation  function  in  time 
domain 
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b.  Response  Evaluation 

-  Create  response  set 

-  initiate  solve 

c.  Save  transient  response  at  a-set  dof  to  a  file  called  'name.unv' 
B.  ACCESS  I-DEAS  TEST  INTERFACE  SOFTWARE 

1.  Correlation  Task 

a.  import  'name.unv'  file  and  change  to  an  Attached  Data  File 
(ADF)  called  'name.ash' 

b.  select  mode  shape  to  work  with 

-  choose  'name.ash' 

-  select  substructure  component 

2.  Post-Processing  Task 

a.  start  animation 
note:  if  using  experimental  data,  once  solve  is  initiated  in  the  Simulation 
Interface  Software,  the  Guyan  transformation  matrix  is  created.  Proceed 
directly  to  Test  Interface  and  load  the  .unv  file  of  a-set  time  histories 


40 


APPENDIX  B.  MATLAB  CODE 


A.        MAIN  FE  PROGRAM 

% 

%         Program  will  solve  for  the  system  mass  and  stiffness 

%         matrices  of  a  plate.  It  uses  four-node  isoparametric 

%         elements  of  the  shear-deformable  displacement 

%         formulation.  The  system  mass  and  stiffness  matrices 

%         are  saved  to  a  file  called  sys_mat.mat 

% 

%  Written  by  Scott  Waltermire 

% 

% 

% 

%         Based  on  EX1071.m 

%         Prof.  Young  Kwon 

% 

% 

% 

% 

% 

% 

%  Variable  descriptions 

%  k  =  element  matrix 

%  kb  =  element  matrix  for  bending  stiffness 

%  ks  =  element  matrix  for  shear  stiffness 

%  f  =  element  vector 

%  kk  =  system  matrix 

%  ff  =  system  vector 

%  disp  =  system  nodal  displacement  vector 

%  gcoord  =  coordinate  values  of  each  node 

%  nodes  =  nodal  connectivity  of  each  element 

%  index  =  a  vector  containing  system  dofs  associated  with  each  element 

%  pointb  =  matrix  containing  sampling  points  for  bending  term 

%  weightb  =  matrix  containing  weighting  coefficients  for  bending  term 

%  points  =  matrix  containing  sampling  points  for  shear  term 

%  weights  =  matrix  containing  weighting  coefficients  for  shear  term 

%  bcdof  =  a  vector  containing  dofs  associated  with  boundary  conditions 

%  bcval  =  a  vector  containing  boundary  condition  values  associated  with 

%  the  dofs  in  'bcdof 

%  kinmtpb  =  matrix  for  kinematic  equation  for  bending 

%  matmtpb  =  matrix  for  material  property  for  bending 

%  kinmtps  =  matrix  for  kinematic  equation  for  shear 

%  matmtps  =  matrix  for  material  property  for  shear 

% 

% 

% 
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% 

%  input  data  for  control  parameters 
% 

clear 


hplate=100; 

lplate=100; 

hele=20; 

lele=20; 

a_ele=hele*lele; 

nel=25; 

nnel=4; 

ndof=3; 

nnode=36; 

sdof=nnode*ndof; 

edof=nnel*ndof; 

emodule=30e6; 

rho=0.284/386.4; 

poisson=0.3; 

t=0.125; 

nglxb=2;  nglyb=2; 

nglb=nglxb*nglyb; 

nglxs=l;  nglys=l; 

ngls=nglxs*nglys; 


%  height  of  plate 

%  length  of  plate 

%  element  height 

%  element  length 

%  element  area 

%  number  of  elements 

%  number  of  nodes  per  element 

%  number  of  dofs  per  node 

%  total  number  of  nodes  in  system 

%  total  system  dofs 

%  degrees  of  freedom  per  element 

%  elastic  modulus 

%  density 

%  Poisson's  ratio 

%  plate  thickness 

%  2x2  Gauss-Legendre  quadrature  for  bending 

%  number  of  sampling  points  per  element  for  bending 

%  lxl  Gauss-Legendre  quadrature  for  shear 

%  number  of  sampling  points  per  element  for  shear 


%- 

% 
% 
% 
% 
% 
% 
% 
%- 


input  data  for  nodal  coordinate  values 
gcoord(i,j)  where  i->node  no.  and  j->x  and  y 
size(gcoord)  =  num_nodes  x  2 


note:    if  the  elements  used  are  not  rectangular 

in  shape,  the  values  of  gcoord  and  nodes 
must  be  manually  entered  by  the  user 


count=0; 

x  =  0:lele:lplate; 

y  =  0:hele:hplate; 

gnode  =  zeros(length(y),length(x)); 

gcoord=zeros(length(x)*length(y),2); 


forn=l:length(y) 

for  j=l:length(x) 

count  =  count+1; 

gcoord(count,l:2)  =  gcoord(count,l:2)  +  [x(j)  y(n)]; 

gnode(nj)  =  gnode(n,j)+count; 
end 
end 
end 
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% 
% 
% 
% 


nodes  is  the  nodal  connectivity  of  the  model 
applied  in  the  counter-clockwise  direction 


nodes=[]; 
count=0; 
forn=l:length(y)-l 

forj=l:length(x)-l 

count  =  count+l; 

node(count,l:4)=[gnode(n,j)  gnode(n,j+l)  gnode(n+l,j+l)  gnode(n+l  j)]; 

nodes  =  [nodes;node(count,l:4)]; 

end 
end 


% 

%  initialization  of  matrices  and  vectors 
% 


ff=zeros(sdof,l); 
kk=zeros(sdof,sdof); 
mm  =  zeros(sdof,sdof); 
disp=zeros(sdof ,  1 ); 
index=zeros(edof,  1); 
kinmtpb=zeros(3  ,edof) ; 

matmtpb=zeros(3,3); 

kinmtps=zeros(2,edof); 

matmtps=zeros(2,2); 


%  system  force  vector 

%  system  k  matrix 

%  system  m  matrix 

%  system  displacement  vector 

%  index  vector 

%  kinematic  matrix  for  bending 

%  constitutive  matrix  for  bending 

%  kinematic  matrix  for  shear 

%  constitutive  matrix  for  shear 


% 

%  computation  of  element  matrices  and  vectors  and  their  assembly 
% 


% 

%  for  bending  stiffness 

% 

[pointb,weightb]=feglqd2(nglxb,nglyb); 
matmtpb=fematiso(  1  ,emodule,poisson)*tA3/12; 


%  sampling  points  &  weights 
%  bending  material  property 


43 


% 

%  for  shear  stiffness 
% 


[points,weights]=feglqd2(nglxs,nglys); 
shearm=0.5*emodule/(  1 .0+poisson); 
shcof=5/6; 
matmtps=shearm*shcof*t*[l  0;  0  1]; 


%  sampling  points  &  weights 
%  shear  modulus 
%  shear  correction  factor 
%  shear  material  property 


for  iel=l:nel 

fori=l:nnel 
nd(i)=nodes(iel,i); 
xcoord(i)=gcoord(nd(i),  1 ); 
ycoord(i)=gcoord(nd(i),2); 
end 

k=zeros(edof ,edof) ; 
m=zeros(edof,edof); 
kb=zeros(edof  ,edof) ; 
ks=zeros(edof  ,edof) ; 


%  loop  for  the  total  number  of  elements 


%  extract  connected  node  for  (iel)-th  element 
%  extract  x  value  of  the  node 
%  extract  y  value  of  the  node 


%  initialize  element  stiffness  matrix  to  zero 
%  initialize  element  mass  matrix  to  zero 
%  initialization  of  bending  matrix  to  zero 
%  initialization  of  shear  matrix  to  zero 


% 

%  numerical  integration  for  bending  term 
% 


for  intx=l:nglxb 
x=pointb(intx,l); 
wtx=weightb(intx,  1 ); 
for  inty=l:nglyb 
y=pointb(inty,2); 
wty=weightb(inty,2) ; 

[shape,dhdr,dhds]=feisoq4(x,y); 


jacob2=fejacob2(nnel,dhdr,dhds,xcoord,y  coord); 

detj  acob=det(j  acob2) ; 

invjacob=inv(jacob2); 


%  sampling  point  in  x-axis 
%  weight  in  x-axis 

%  sampling  point  in  y-axis 
%  weight  in  y-axis 

%  compute  shape  functions  and 
%  derivatives  at  sampling  point 

%  compute  Jacobian 

%  determinant  of  Jacobian 

%  inverse  of  Jacobian  matrix 


[dhdx,dhdy]=federiv2(nnel,dhdr,dhds,invjacob);       %  derivatives  w.r.t. 

%  physical  coordinate 


kinmtpb=fekinepb(nnel,dhdx,dhdy); 


%  bending  kinematic  matrix 


% 

%  compute  bending  and  mass  element  matrix 
% 
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kb=kb+kinmtpb'*matmtpb*kinmtpb*wtx*wty*detjacob; 


% 

% 
% 

% 

Q  = 


Create  a  matrix  of  shape  functions:  Q 

Ensure  each  row  vector  in  the  matrix  is  properly  weighted:  H(i) 


[shape(l)  0  0  shape(2)  0  0  shape(3)  0  0  shape(4)  0  0;... 
0  shape(l)  0  0  shape(2)  0  0  shape(3)  0  0  shape(4)  0; 
0  0  shape(l)  0  0  shape(2)  0  0  shape(3)  0  0  shape(4)]; 

Hl=shape(l)*Q; 

HI  =  [l/12*rho*tA3*Hl(l,:);l/12*rho*tA3*Hl(2,:);rho*t*Hl(3,:)J; 

H2  =  shape(2)*Q; 

H2  =  [l/12*rho*tA3*H2(l,:);l/12*rho*tA3*H2(2,:);rho*t*H2(3,:)]; 

H3  =  shape(3)*Q; 

H3  =  [l/12*rho*tA3*H3(l,:);l/12*rho*tA3*H3(2,:);rho*t*H3(3,:)]; 

H4  =  shape(4)*Q; 

H4  =  [l/12*rho*tA3*H4(l,:);l/12*rho*tA3*H4(2,:);rho*t*H4(3,:)]; 


% 

H  = 

% 
% 
% 


Create  the  elemental  mass  matrix  in  terms  of  shape  functions 
[H1;H2;H3;H4]; 

Solve  for  the  elemental  mass  matrix 


m  =  m+H*wtx*wty*detjacob; 


end 
end 


%  end  of  numerical  integration  loop  for  bending  term 

%  and  mass  term 


% 

%  numerical  integration  for  shear  term 
% „__ 

for  intx=l:nglxs 
x=points(intx,l); 
wtx=weights(intx,  1 ); 
for  inty=l:nglys 
y=points(inty,2); 
wty=weights(inty,2) ; 

[shape,dhdr,dhds]=feisoq4(x,y); 


jacob2=fejacob2(nnel,dhdr,dhds,xcoord,ycoord); 
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%  sampling  point  in  x-axis 
%  weight  in  x-axis 

%  sampling  point  in  y-axis 
%  weight  in  y-axis 

%  compute  shape  functions  and 
%  derivatives  at  sampling  point 

%  compute  Jacobian 


detj acob=det(j acob2);  %  determinant  of  Jacobian 

invjacob=inv(jacob2);  %  inverse  of  Jacobian  matrix 

[dhdx,dhdy]=federiv2(nnel,dhdr,dhds,invjacob);       %  derivatives  w.r.t. 

%  physical  coordinate 

kinmtps=fekineps(nnel,dhdx,dhdy,shape);  %  shear  kinematic  matrix 

% 

%  compute  shear  element  matrix 
% 

ks=ks+kinmtps'*matmtps*kinmtps*wtx*wty*detjacob; 

end 

end  %  end  of  numerical  integration  loop  for  shear  term 

% 

%  compute  element  matrix 
% 

k=kb+ks; 

index=feeldof(nd,nnel,ndof);  %  extract  system  dofs  associated  with  element 

kk=feasmbl  1  (kk,k,index);  %  assemble  element  matrices  into  system 

mm=f easmbl  1  (mm,m,index) ;  %  matrices 

end 

%  Apply  Springs  in  z-direction 

kk(3,3)  =  kk(3,3)  +  100; 
kk(18,18)  =  kk(18,18)+  100; 
kk(93,93)  =  kk(93,93)  +  100; 
kk(108,108)  =  kk(108,108)  +  100; 

save  sys_mat  mm  kk 
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B. 


FUNCTIONS  CALLED  BY  MAIN  FE  PROGRAM 


function  [point2,weight2]=feglqd2(nglx,ngly) 

% 

%  Purpose: 

%     determine  the  integration  points  and  weighting  coefficients 

%      of  Gauss-Legendre  quadrature  for  two-dimensional  integration 

% 

%  Synopsis: 

%      [point2,weight2]=feglqd2(nglx,ngly) 

% 

%  Variable  Description: 

%      nglx  -  number  of  integration  points  in  the  x-axis 

%      ngly  -  number  of  integration  points  in  the  y-axis 

%     point2  -  vector  containing  integration  points 

%     weight2  -  vector  containing  weighting  coefficients 

% 

%  determine  the  largest  one  between  nglx  and  ngly 

if  nglx  >  ngly 

ngl=nglx; 
else 

ngl=ngly; 
end 

%  initialization 

point2=zeros(ngl,2) ; 
weight2=zeros(ngl,2); 

%  find  corresponding  integration  points  and  weights 


[pointx,  weightx]=feglqd  1  (nglx) ; 
[pointy,weighty]=feglqdl(ngly); 

%  quadrature  for  two-dimension 


%  quadrature  rule  for  x-axis 
%  quadrature  rule  for  y-axis 


for  intx=l:nglx 

point2(intx,  1  )=pointx(intx); 

weight2(intx,  1  )=weightx(intx); 
end 


%  quadrature  in  x-axis 


for  inty=l:ngly 

point2(inty,2)=pointy(inty); 

weight2(inty,2)=weighty(inty); 
end 


%  quadrature  in  y-axis 
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function  [matmtrx]=fematiso(iopt,elastic,poisson) 

% 

%  Purpose: 

%     determine  the  constitutive  equation  for  isotropic  material 

% 

%  Synopsis: 

%      [matmtrx]=fematiso(iopt,elastic,poisson) 

% 

%  Variable  Description: 

%      elastic  -  elastic  modulus 

%      poisson  -  Poisson's  ratio 

%      iopt=l  -  plane  stress  analysis 

%      iopt=2  -  plane  strain  analysis 

%      iopt=3  -  axisymmetric  analysis 

%      iopt=4  -  three  dimensional  analysis 

% 


'o- 


if  iopt==l         %  plane  stress 
matmtrx=  elastic/(l-poisson*poisson)*  ... 
[1  poisson  0;  ... 
poisson   1   0;  ... 
0  0  (l-poisson)/2]; 

elseif  iopt=2      %  plane  strain 
matmtrx=  elastic/((  l+poisson)*(  1  -2*poisson))* 
[(1 -poisson)  poisson  0; 
poisson  (1 -poisson)  0; 
0  0  (l-2*poisson)/2]; 

elseif  iopt==3      %  axisymmetry 
matmtrx=  elastic/((  1  +poisson)  *  ( 1  -2*poisson))  * 
[(1 -poisson)  poisson  poisson  0; 
poisson  (1 -poisson)    poisson  0; 
poisson  poisson  (1 -poisson)    0; 
0     0     0    (l-2*poisson)/2]; 

else      %  three-dimension 
matmtrx=  elastic/((l+poisson)*(l-2*poisson))* 
[(1 -poisson)  poisson  poisson    0    0     0; 
poisson  (1 -poisson)    poisson    0    0     0; 
poisson  poisson  (1 -poisson)     0    0     0; 
0     0     0     (l-2*poisson)/2    0     0; 
0     0     0     0     (l-2*poisson)/2   0; 
0     0    0     0     0    (l-2*poisson)/2]; 

end 
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function  [shapeq4,dhdrq4,dhdsq4]=feisoq4(rvalue,svalue) 

% 

%  Purpose: 

%     compute  isoparametric  four-node  quadilateral  shape  functions 

%     and  their  derivatves  at  the  selected  (integration)  point 

%     in  terms  of  the  natural  coordinate 

% 

%  Synopsis: 

%      [shapeq4,dhdrq4,dhdsq4]=feisoq4(rvalue,svalue) 

% 

%  Variable  Description: 

%      shapeq4  -  shape  functions  for  four-node  element 

%      dhdrq4  -  derivatives  of  the  shape  functions  w.r.t.  r 

%      dhdsq4  -  derivatives  of  the  shape  functions  w.r.t.  s 

%      rvalue  -  r  coordinate  value  of  the  selected  point 

%      svalue  -  s  coordinate  value  of  the  selected  point 

% 

%  Notes: 

%      1st  node  at  (-1,-1),  2nd  node  at  (1,-1) 

%      3rd  node  at  (1,1),  4th  node  at  (-1,1) 

% 

%  shape  functions 

shapeq4(l)=0.25*(l-rvalue)*(l-svalue); 
shapeq4(2)=0.25*(l+rvalue)*(l  -svalue); 
shapeq4(3)=0.25*(l+rvalue)*(l+svalue); 
shapeq4(4)=0. 25  *  ( 1  -rvalue)  *  ( 1  +s  value) ; 

%  derivatives 

dhdrq4(l)=-0.25*(l-svalue); 
dhdrq4(2)=0.25*(l-svalue); 
dhdrq4(3)=0.25*(l+svalue); 
dhdrq4(4)=-0.25*(l+svalue); 

dhdsq4(l)=-0.25*(l-rvalue); 
dhdsq4(2)=-0.25*(l+rvalue); 
dhdsq4(3)=0.25*(l+rvalue); 
dhdsq4(4)=0.25*  ( 1  -rvalue); 


49 


function  (jacob2]=fejacob2(nnel,dhdr,dhds,xcoord,ycoord) 

% 

%  Purpose: 

%     determine  the  Jacobian  for  two-dimensional  mapping 

% 

%  Synopsis: 

%      [jacob2]=fejacob2(nnel,dhdr,dhds,xcoord,ycoord) 

% 

%  Variable  Description: 

%  jacob2  -  Jacobian  for  one-dimension 

%  nnel  -  number  of  nodes  per  element 

%  dhdr  -  derivative  of  shape  functions  w.r.t.  natural  coordinate  r 

%  dhds  -  derivative  of  shape  functions  w.r.t.  natural  coordinate  s 

%  xcoord  -  x  axis  coordinate  values  of  nodes 

%  ycoord  -  y  axis  coordinate  values  of  nodes 

% 

jacob2=zeros(2,2); 

for  i=l:  nnel 

jacob2(  1 , 1  )=jacob2(  1 , 1  )+dhdr(i)*xcoord(i) ; 

jacob2(l,2)=jacob2(l,2)+dhdr(i)*ycoord(i); 

jacob2(2,l)=jacob2(2,l)+dhds(i)*xcoord(i); 

jacob2(2,2)=jacob2(2,2)+dhds(i)*ycoord(i); 

end 


function  [dhdx,dhdy]=federiv2(nnel,dhdr,dnds,invjacob) 

% 

%  Purpose: 

%  determine  derivatives  of  2-D  isoparametric  shape  functions  with 

%  respect  to  physical  coordinate  system 

% 

%  Synopsis: 

%  [dhdx,dhdy]=federiv2(nnel,dhdr,dnds,invjacob) 

% 

%  Variable  Description: 

%      dhdx  -  derivative  of  shape  function  w.r.t.  physical  coordinate  x 

%      dhdy  -  derivative  of  shape  function  w.r.t.  physical  coordinate  y 

%     nnel  -  number  of  nodes  per  element 

%     dhdr  -  derivative  of  shape  functions  w.r.t.  natural  coordinate  r 

%      dhds  -  derivative  of  shape  functions  w.r.t.  natural  coordinate  s 

%     invjacob  -  inverse  of  2-D  Jacobian  matrix 

% 


for  i=l:nnel 

dhdx(i)=invjacob(  1 , 1  )*dhdr(i)+invjacob(  1 ,2)*dhds(i); 

dhdy(i)=invjacob(2, 1  )*dhdr(i)+invjacob(2,2)*dhds(i); 

end 
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function  [kinmtpb]=fekinepb(nnel,dhdx,dhdy) 

% 

%  Purpose: 

%  determine  the  kinematic  matrix  expression  relating  bending  curvatures 

%  to  rotations  and  displacements  for  shear  deformable  plate  bending 

% 

%  Synopsis: 

%  [kinmtpb]=fekinepb(nnel,dhdx,dndy) 

% 

%  Variable  Description: 

%  nnel  -  number  of  nodes  per  element 

%  dhdx  -  derivatives  of  shape  functions  with  respect  to  x 

%  dhdy  -  derivatives  of  shape  functions  with  respect  to  y 

% 

fori=l:nnel 

il=(i-l)*3+l; 

i2=il+l; 

i3=i2+l; 

kinmtpb(  1  ,i  1  )=dhdx(i); 

kinmtpb(2,i2)=dhdy(i); 

kinmtpb(3  ,i  1  )=dhdy  (i) ; 

kinmtpb(3,i2)=dhdx(i); 

kinmtpb(3,i3)=0; 

end 


function  [kk]=feasmbll(kk,k,index) 

% 

%  Purpose: 

%  Assembly  of  element  matrices  into  the  system  matrix 

% 

%  Synopsis: 

%  [kk]=feasmbll(kk,k,index) 

% 

%  Variable  Description: 

%  kk  -  system  matrix 

%  k  -  element  matri 

%  index  -  d.o.f.  vector  associated  with  an  element 

% 


edof  =  length(index); 
for  i=l:edof 
ii=index(i); 
for  j=l:edof 
jj=index(j); 

kk(ii,jj)=kk(ii,jj)+k(i,j); 
end 
end 
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function  [kinmtps]=fekineps(nnel,dhdx,dhdy,shape) 

% 

%  Purpose: 

%      determine  the  kinematic  matrix  expression  relating  shear  strains 

%      to  rotations  and  displacements  for  shear  deformable  plate  bending 

% 

%  Synopsis: 

%      [kinmtps]=fekineps(nnel,dhdx,dhdy,shape) 

% 

%  Variable  Description: 

%     nnel  -  number  of  nodes  per  element 

%      dhdx  -  derivatives  of  shape  functions  with  respect  to  x 

%      dhdy  -  derivatives  of  shape  functions  with  respect  to  y 

%      shape  -  shape  function 

% 

for  i=l:nnel 

il=(i-l)*3+l; 

i2=il+l; 

i3=i2+l; 

kinmtps(  1  ,i  1  )=-shape(i); 

kinmtps(  1  ,i3)=dhdx(i); 

kinmtps(2,i2)=-shape(i); 

kinmtps(2,i3)=dhdy(i); 

end 

function  [index]=feeldof(nd,nnel,ndof) 

% 

%  Purpose: 

%      Compute  system  dofs  associated  with  each  element 

% 

%  Synopsis: 

%      [index] =feeldof(nd,nnel,ndof) 

% 

%  Variable  Description: 

%     index  -  system  dof  vector  associated  with  element  "iel" 

%     nd  -  nodal  vector  whose  system  dofs  are  to  be  determined 

%      nnel  -  number  of  nodes  per  element 

%      ndof  -  number  of  dofs  per  node 

% 

edof  =  nnel*ndof; 
k=0; 

for  i=l:  nnel 
start  =  (nd(i)-l)*ndof; 
forj=l:ndof 
k=k+l; 

index(k)=start+j; 
end 
end 
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MAIN  POST-PROCESSING  PROGRAM 


clear 


% 

% 

%  This  program  will  load  the  system  mass  and  stiffness  matrices 

%  from  file  sys_mat.mat.  It  solves  the  eigenvalue  problem 

%  to  get  the  natural  frequencies  and  the  modeshapes. 

%  Convolution  is  used  to  solve  for  the  modal  solution  q(t). 

%  It  solves  for  x(t)  at  all  nodes  then  calls  set.mat  to 

%  pick  out  the  response  at  the  analysis  set.  From  this  it 

%  back-expands  the  aset  into  the  full  solution  using  Guyan  and 

%  IRS  transformation  matrices.  It  can  plot  the  modeshapes, 

%  time-histories  of  each  solution  and  the  error  between 

%  the  actual  response  and  the  Guyan/IRS  models. 

% 

%  Written  by  Scott  Waltermire 

%  Get  the  global  mass  and  stiffness  matrices 

load  sys_mat; 
sdof=108; 

%  Create  Blast  Forcing  Function 

plotit  =  1;  %  plotit=l  plots  forcing  function 

Fo  =  2000;  %  amplitude  of  forcing  function 

time  =  0:.01:3;  %  time  step 

[f_of_t]  =  fBlastForcing(Fo,time,plotit);  %  function  to  solve  for  f(t) 

F  =  zeros(sdof,length(time));  %  initialize  system  force  vector 

F(63,:)  =  f_of_t;  %  place  f(t)  at  the  proper  node 


%  Solve  and  Sort  Eigenvalues/Modeshapes 

[phi,lam]=eig(mm\kk) ; 

mtilda=phi'*mm*phi;  %  mass  normalize 

%  loop  to  mass  normalize  the  modeshapes 

for  i=l:length(mm) 
phi(:,i)=phi(:,i)*l/(sqrt(mtilda(i,i))); 
end 

%  loop  to  pick  out  the  natural  frequencies 

for  j= 1  :length(mm) 
ev(j)=lam(j,j); 
end 

[lam,p]=sort(ev);  %  sort  the  natural  frequencies 
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phistar=phi; 

%  loop  to  sort  the  eigenvectors  in  the  same  manner 
%  as  the  natural  frequencies 

for  k=l  ;length(mm) 
phi(:,k)=phistar(:,p(k)); 
end 

wn=real(sqrt(lam)); 

omega=real((sqrt(lam)/(2*pi)))';  %  convert  wn  to  hertz 

phi  =  real(phi); 

%  create  modal_F  and  use  proportional  damping 

modal_F  =  phi'*F; 
zeta  =  .02; 

%  Convolution  to  solve  for  q(t)  using  Trapezoidal  rule 

dt  =  .001; 
fori  =  1:108 

wd(i)  =  wn(i)*sqrt(l-zetaA2); 
h  =  l/wd(i)*exp(-zeta*wn(i)*time).*sin(wd(i)*time); 
F  =  modal_F(i,:); 
[convXY]  =  fconvTrap(F,h,dt); 
q(i,:)  =  convXY; 
end 

%  Plot  Modeshapes 

num_mode_shape=  10; 
for  i=l  :num_mode_shape 

phi_mode=real(phi(:  ,i)) ; 

phil=(phi_mode(3:3:18))'; 

phi2=(phi_mode(21 :3:36))'; 

phi3=(phi_mode(39:3:54))'; 

phi4=(phi_mode(57:3:72))'; 

phi5=(phi_mode(75:3:90))'; 

phi6=(phi_mode(93:3: 108))'; 

phi_mesh=[phi  1  ;phi2;phi3;phi4;phi5;phi6]; 

surfc(phi_mesh) 

shading  interp 

grid 

pause 
end 
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%  Convert  back  to  physical  Coordinates 

resp  =  phi*q; 

%  Solve  for  T_static  and  T_irs 

load  set;  %  aset/oset  file 

%         functions  to  solve  for  Guyan  and  IRS 
%         transformation  matrices. 

[kstat,mstat,T_static]  =  fstatic_tam(kk,mm,oset,aset); 
[kirs,mirs,T_irs]  =  firs_tam(kk,mm,oset,aset); 


% 

%  Load  Guyan  Transformation  matrix 

%  and  determine  xa  and  xo 

% 

for  i  =  l:length(aset) 

xa_g(i,:)  =  resp(aset(i),:); 
end 

%  Back  expand  the  guyan  aset  time  histories 
%  into  the  full  response 

x_g  =  T_static*xa_g; 

xo_g  =  x_g(l:sdof-length(aset),:); 

xa_g  =  x_g(sdof-length(aset)+l:sdof,:); 

%  Now  re-sort  xa_g  and  xo_g 

fori=l:length(aset) 

xol_g  =  xo_g(l:aset(i)-l,:); 

xo2_g  =  xo_g(aset(i):length(xo_g(:,l)),:); 

xol_g  =  [xol_g;xa_g(i,:)]; 

xo_g  =  [xol_g;xo2_g]; 
end 
x_g  =  xo_g; 


% 

%  Load  IRS  Transformation  matrix 

%  and  determine  xa  and  xo 

% 


fori  =  l:length(aset) 

xa_irs(i,:)  =  resp(aset(i),:); 
end 
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%  Back  expand  the  irs  aset  time  histories 
%  into  the  full  response 

x_irs  =  T_irs*xa_irs; 

xo_irs  =  x_irs(l:sdof-length(aset),:); 

xa_irs  =  x_irs(sdof-length(aset)+l:sdof,:); 

%  Now  re-sort  xa_irs  and  xo_irs 

fori=l:length(aset) 

xol_irs  =  xo_irs(l:aset(i)-l,:); 

xo2_irs  =  xo_irs(aset(i):length(xo_irs(:,l)),:); 

xol_irs  =  [xol_irs;xa_irs(i,:)]; 

xo_irs  =  [xol_irs;xo2_irs]; 
end 
x_irs  =  xo_irs; 


%  Ignore  all  rotations  and  plot  response 
%  in  the  z-direction 

resp  =  resp(3:3:sdof,:); 
resp_g  =  x_g(3:3:sdof,:); 
resp_irs  =  x_irs(3:3:sdof,:); 

for  i=  1:36 

figure(l) 
plot(time,resp(i,:),'r-.',time,resp_g(i,:),'y'); 
xlabel('Time  (sec)'); 
ylabel('Displacement  (in)'); 
title('Time  History'); 

figure(2) 
plot(time,resp(i,:),'r-.',time,resp_irs(i,:),'y'); 
xlabel('Time  (sec)'); 
ylabel('Displacement  (in)'); 
title('Time  History'); 

pause 
end 

close(l) 
close(2) 
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% 

%         Find  the  error  between  the  true  response  and 
%         the  Guyan/IRS  back-expanded  response 
% 

for  i  =  1  :length(time) 

xg_err(:,i)  =  abs(resp_g(:,i)  -  resp(:,i)); 

xirs_err(:,i)  =  abs(resp_irs(:,i)  -  resp(:,i)); 
end 
fori  =1:36 

xg_error(i)  =  max(xg_err(i,:)); 

xirs_error(i)  =  max(xirs_err(i,:)); 
end 

xg_errorl  =  xg_error(l:6); 

xg_error2  =  xg_error(7:12); 

xg_error3  =  xg_error(  13:18); 

xg_error4  =  xg_error(  19:24); 

xg_error5  =  xg_error(25:30); 

xg_error6  =  xg_error(31:36); 

xg_error  =  [xg_errorl  ;xg_error2;xg_error3;xg_error4;xg_error5;xg_error6]; 

figure(l) 

surf(xg_error); 

shading  interp 

grid 

xirs_errorl  =  xirs_error(l:6); 

xirs_error2  =  xirs_error(7:12); 

xirs_error3  =  xirs_error(  13:18); 

xirs_error4  =  xirs_error(  19:24); 

xirs_error5  =  xirs_error(25:30); 

xirs_error6  =  xirs_error(31:36); 

xirs_error  =  [xirs_erTorl;xirs_error2;xirs_error3;... 

xirs_error4;xirs_error5  ;xirs_error6] ; 
figure(2) 
surf(xirs_error); 
shading  interp 
grid 


save  modal_info  resp  resp_g  resp_irs  time 
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D.         FUNCTIONS  CALLED  BY  MAIN  PROGRAM 

function  [f_of_t]  =  fBlastForcing(Fo,time,plotit); 

% 

%     This  function  returns  a  forcing  function  which  is 

%     a  "blast"  function. 

% 

%  F(t)  =  Fo  *  (  exp(-at)  -  exp(-bt) ) 

% 

%     where  a  and  b  are  constants  which  shape  the  blast, 

%     and  Fo  is  the  amplitude  of  the  blast. 

% 

%     The  variable  "plotit"  is  a  switch  which  if  set  =  1  will 

%     cause  the  f(t)  to  be  plotted,  if  set  to  anything  else 

%     will  not  plot. 

%  written  by  J.H.  Gordis 

% 

%  Choices:  sine    blst     step 

type  =  'blst'; 


% 

if  type  ==  'blst'; 

disp('  Blast  forcing  used...') 

a=  100.0; 

b  =  300.0; 

f_of_t  =  Fo  *  (  exp(-a*time)  -  exp(-b*time) ); 

elseif  type  ==  'step'; 

%        dispC  Step  forcing  used...') 
f_of_t  =  Fo  *  ones(size(time)); 

elseif  type  ==  'sine'; 

%        dispC  Sine  forcing  used...') 
W  =  5;  %  Hertz 
f_of_t  =  Fo  *  sin(2*pi*W*time); 

end; 

if  plotit  ==  1 ; 

figure(9) 

plot(time,f_of_t)  ;grid 

pause 

elf 
end 

%  End  function. 
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function  [convXY]  =  fConvTrap(x,y,dt); 

% 

%  Usage:    [convXY]  =  fConvTrap(x,y,dt); 

% 

%  This  function  performs  the  convolution  integral  calculation 

% 

%  tau=t 

%      convXY(t)  =     S      x(t-tau)  y(tau)  d  tau 

%  tau=0 

% 

%  using  the  trapezoidal  rule. 

%  The  function  is  passed  the  vectors  x(t)  and  y(t) 

%  and  the  sample  spacing  (time  step)  dt. 

% 

%     x       =  vector  of  length  n 

%     y       =  vector  of  length  n 

%     dt     =  sample  spacing  (i.e.  delta_t). 

%  written  by  J.H.  Gordis 

% 

if  size(x)  ==  size(y);         %  Vectors  the  same  size.  Perform  convolution. 

% 

convXY  =  zeros(size(x),l); 

for  icnt_step  =  2  :  length(x); 
[wts]  =  fTrapzWts(icnt_step); 
if  size(x,l)  ~=  size(wts,l); 

wts  =  wts'; 
end 
convXY(icnt_step)  =  dt  *  sum(wts  .*  fliplr(x(l:icnt_step))  .*  y(l:icnt_step) ); 
end; 

elseif  size(x)  ~=  size(y)  &  length(x)  ==  length(y);    %  Vectors  are  transposed.  Don't 
perform  convolution. 

% 


disp('  Error  in  fConvTrap.  Transpose  one  vector  and  try  again.') 

disp(sprintf('  Size(x)  %3i',size(x))) 

disp(sprintf('  Size(y)  %3i',size(y))) 

elseif  size(x)  ~=  size(y)  &  length(x)  ~=  length(y);    %  Vectors  different  sizes.  Don't 
perform  convolution. 

% 


disp('  Error  in  fConvTrap.  Vector  not  the  same  size.') 

disp(sprintf('  Size(x)  %3i',size(x))) 

disp(sprintf('  Size(y)  %3i',size(y))) 

end; 
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% 

function  [kstat,mstat,T_static]=fstatic_tam(k,m,oset,aset) 

% 

%  this  function  returns  the  statically  reduced  stiffness 

%  and  mass  matrices,  given  the  unreduced  couterparts. 

%  Care  must  be  taken  that  the  aset  and  oset  vectors  correspond 

%  with  the  existing  arrangement  of  k  and  m. 

%  k  and  m  are  UNPARTTTIONED  matrices. 

% 

%  written  by  J.H.  Gordis 

aset_size=length(aset) ; 

% 

kaa=k(aset,aset); 

kao=k(aset,oset); 

koo=k(oset,oset); 

koa=kao'; 

clear  k; 

k=[koo,koa;kao,kaa]; 

% 

maa=m(aset,aset); 

mao=m(aset,oset) ; 

moo=m(oset,oset) ; 

moa=mao'; 

clear  m; 

m=[moo,moa;mao,maa]; 

% 

t_static=-koo\koa; 

T_static  =  [t_static;eye(aset_size)]; 

% 

kstat=T_static'*k*T_static; 

mstat=T_static'*m*T_static; 

% 

%  end  function  fstatic  tarn 
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% 

function  [kirs,mirs,T_irs]=firs_tam(k,m,oset,aset) 

% 

%  this  function  returns  the  IRS  reduced  stiffness 

%  and  mass  matrices,  given  the  unreduced  couterparts. 

%  Care  must  be  taken  that  the  aset  and  oset  vectors  correspond 

%  with  the  existing  arrangement  of  k  and  m. 

%  k  and  m  are  UNPARTTTIONED  matrices. 

%  written  by  J.H.  Gordis 

aset_size=length(aset) ; 

% 

kaa=k(aset,aset); 

kao=k(aset,oset); 

koo=k(oset,oset); 

koa=kao'; 

clear  k; 

k=[koo,koa;kao,kaa]; 

% 

maa=m(aset,aset) ; 

mao=m(aset,oset) ; 

moo=m(oset,oset) ; 

moa=mao'; 

clear  m; 

m=[moo,moa;mao,maa] ; 

% 

t_static=-koo\koa; 

T_static  =  [t_static;eye(aset_size)]; 

% 

kstat=T_static'*k*T_static; 

mstat=T_static'*m*T_static; 

% 

tirs=t_static+inv(koo)*(moa+moo*t_static)*inv(mstat)*kstat; 

T_irs= [tirs  ;ey  e(aset_size)] ; 

% 

kirs=T_irs'*k*T_irs; 

mirs=T_irs'*m*T_irs; 

% 

%  end  function  firs  tarn 
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% 

% 

%  File  holding  the  dof  for  the  analysis  set  (aset) 

%  and  the  omitted  set  (oset) 

% 

%  Saved  to  a  file  called  set.mat 

% 

% 

aset  =[3  18  93  108]; 

oset=[l  245  67  8  9  10  11  12  13  14  15  16  17... 

19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  ... 

40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61 

62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80  81  82  ... 

83  84  85  86  87  88  89  90  91  92  94  95  96  97  98  99  100  101  102  ... 

103  104  105  106  107]; 

save  set  aset  oset 


%         Program  will  animate  the  full-order  response 
%         by  loading  a  .mat  file  called  modal-info 

clear; 

load  modal_info; 

clear  resp_g  resp_irs 

%  make  movie 

movie_fig  =  figure('position',[100  200  300  200]); 
M  =  moviein(200); 
[x,y]  =  meshgrid([-5:2:5]); 
fori  =1:200 

resp_t  =  resp(:,i); 

respl  =  (resp_t(l:6))'; 

resp2  =  (resp_t(7:12))'; 

resp3  =  (resp_t(13:18))'; 

resp4  =  (resp_t(  19:24))'; 

resp5  =  (resp_t(25:30))'; 

resp6  =  (resp_t(31:36))'; 

z  =  [respl  ;resp2;resp3;resp4;resp5;resp6]; 
surf(x,y,z); 
shading  interp 
grid; 

axis([-5  5-5  5  -.4  .6]); 
view([45  45  10]) 
zlabel('disp  (in.)'); 
M(:,i)  =  getframe; 
end 
pause 
movie(M,0); 
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%         Program  will  animate  the  plate  response  found 
%         by  the  static  reduction  method.  It  loads  the  .mat 
%         file  called  modal_info 

clear; 

load  modal_info; 

clear  resp  resp_irs 
resp  =  resp_g; 

%  make  movie  of  guyan  response 

movie_fig  =  figure('position',[100  200  300  200]); 

M  =  moviein(200); 

[x,y]  =  meshgrid([-5:2:5]); 

count  =  1 ; 
fori  =1:200 

resp_t  =  resp(:,i); 

respl  =  (resp_t(l:6))'; 

resp2  =  (resp_t(7:12))'; 

resp3  =  (resp_t(13:18))'; 

resp4  =  (resp_t(  19:24))'; 

resp5  =  (resp_t(25:30))'; 

resp6  =  (resp_t(31:36))'; 

z  =  [respl  ;resp2;resp3;resp4;resp5;resp6]; 
surf(x,y,z); 
shading  interp 
grid; 

axis([-5  5-5  5  -.4  .6]); 
view([45  45  10]) 
zlabel('disp  (in.)'); 

title(' Animation  by  Guyan  Back  Expansion'); 
M(:, count)  =  getframe; 
count  =  count+l; 
end 

pause 

movie(M,0); 
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% 

%  Program  to  animate  the  plate  response  found 

%  by  the  IRS  reduction  method.  Loads  the  .mat  file 

%  modal_info 

% 

clear; 

load  modal_info; 

clear  resp  resp_g 
resp  =  resp_irs; 

%  make  movie  of  guyan  response 

movie_fig  =  figure(*position',[100  200  300  200]); 

M  =  moviein(200); 

[x,y]  =  meshgrid([-5:2:5]); 

count  =  1 ; 
fori  =1:200 

resp_t  =  resp(:,i); 

respl  =  (resp_t(l:6))'; 

resp2  =  (resp_t(7:12))'; 

resp3  =  (resp_t(13:18))'; 

resp4  =  (resp_t(  19:24))'; 

resp5  =  (resp_t(25:30))' 

resp6  =  (resp_t(31:36))' 

z  =  [respl  ;resp2;resp3;resp4;resp5;resp6]; 
surf(x,y,z); 
shading  interp 
grid; 

axis([-5  5-5  5  -.4  .6]); 
view([45  45  10]) 
zlabel('disp  (in.)'); 

title('Animation  by  IRS  Back  Expansion'); 
M(:, count)  =  getframe; 
count  =  count+1; 
end 

pause 

movie(M,0); 
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